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INTRODUCTION 


The  present  study  considers  both  laminar  and  turbulent  flow  in  curved  ducts, 
pipes,  and  channels  of  constant  cross  sectional  area  and  shape.  The  particular  flows 
considered  here,  as  well  as  similar  flows  in  related  geometries,  are  very  common  in 
internal  flow  applications.  They  are  of  interest  in  connection  with  flow  degradation 
and  their  strong  influence  on  flow  losses  and  heat  transfer  levels. 

The  general  character  of  flow  in  curved  ducts  and  pipes  is  known  to  differ 
fundamentally  from  that  in  straight  flow  geometries,  due  to  the  presence  of  large 
secondary  flows  which  distort  the  primary  flow.  Since  strong  deflections  may  occur 
over  a  short  distance,  such  flows  are  usually  of  a  transition  type  and  seldom  become 
fully  developed  or  assume  any  convenient  similarity  form.  By  analogy  with  external 
flows,  the  flow  often  behaves  as  an  inviscid  flow  in  a- central  core  region,  with 
viscous  effects  limited  to  regions  near  solid  boundaries,  Unlike  most  external  flow 
problems,  however,  the  inviscid  core  region  is  often  not  an  irrotational  potential 
flow  but  is  rotational  with  interaction  between  the  viscous  and  rotational  Inviscid 
flow  regions.  Furthermore,  as  the  flow  passes  through  successive  flow  passages,  new 
viscous  and  thermal  boundary  layers  develop  beneath  previous  boundary  layers,  and  the 
distinction  between  rotational  inviscid  and  viscous  boundary  layer  regions  can  become 
tenuous . 


The  underlying  physical  mechanisms  present  in  flows  of  this  type  are  clearly  j 

elucidated  by  secondary  flow  theory  (reviewed  for  example  by  Horlock  &  Lakshminarayana  [1]).| 
In  its  most  common  and  perhaps  simplest  form,  the  secondary  flow  is  generated  by  turning  ] 
a  primary  flow  in  which  viscous  or  other  forces  upstream  have  produced  a  non-zero  1 

velocity  gradient  normal  to  the  plane  of  curvature.  Fluid  with  above  (/below)  average  j 
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momentum  migrates  to  the  outside  (/inside)  of  the  bend  as  a  result  of  the  radial 
pressure  gradients  produced  by  turning  the  flow.  This  phenomenon  is  quantified  by 
secondary  flow  theory  as  the  generation  of  streamwise  vorticity  from  transverse 
vorticity  which  has  been  produced  upstream.  Although  the  secondary  flow  is  generated 

i 

by  an  inviscid  mechanism,  its  strength  and  subsequent  development  are  influenced  in 
varying  degrees  by  viscous  effects.  In  any  event,  the  secondary  flows  are  often 
quite  large  and  the  flow  patterns  are  thus  complex  and  highly  three-dimensional. 
Two-dimensional  flows  in  curved  channels  do  not  behave  in  this  manner  (provided  the 
flow  remains  two-dimensional),  since  there  can  be  no  streamwise  vorticity  component 
in  a  two-dimensional  flow. 


Previous  Work 


Much  of  the  early  work  on  flow  in  curved  ducts  and  pipe  bends  can  be  traced 
from  the  reviews  of  Hawthrone  [2,  3].  Of  primary  interest  here  are  methods  for 
computing  such  flows  and  experimental  work  which  may  be  useful  in  evaluating  the 
numerical  predictions.  Only  three-dimensional  developing  flows  in  curved  ducts  and 
pipe  bends  are  considered.  Two-dimensional  flows  and  flows  which  are  fully  developed 
and  thus  do  not  vary  with  an  axial  coordinate  are  excluded  from  consideration. 

Flow  in  Curved  Ducts 

Pratap  and  Spalding  [4]  have  considered  turbulent  flow  in  a  strongly  curved  duct 
{  using  their  "partially  parabolic"  calculation  procedure  and  a  two-equation/wall-function 
turbulence  model.  Ghia  and  Sokhey  [5]  have  computed  laminar  flow  in  ducts  of  strong 
curvature  using  a  parabolized  form  of  the  Navier-Stokes  equations.  Kreskovsky,  Briley, 
and  McDonald  [6]  have  recently  applied  an  approximate  initial-value  analysis  for 
viscous  primary  and  secondary  flows  to  laminar  and  turbulent  flow  in  strongly  curved 
ducts,  using  a  one-equation  turbulence  model  with  viscous  sublayer  resolution. 

Humphrey,  Taylor  and  Whitelaw  [7]  have  obtained  laser-Doppler  anemometry  measurements 
of  laminar  flow  in  a  square  duct  of  strong  curvature,  for  the  case  of  fully  developed 
flow  in  the  straight  section  upstream  of  the  bend.  For  comparison,  then  also  performed 
numerical  calculations  for  this  flow  using  a  version  of  the  fully-elliptic  calculation 
procedure  developed  at  Imperial  College  by  Gosman,  Pun,  Patankar  and  Spalding. 

Further  extensive  calculations  including  heat  transfer  effects  have  been  made  recently 
by  Yee  and  Humphrey  [8].  Finally,  Taylor,  Whitelaw  and  Yianneskis  [9]  have  recently 
made  extensive  measurements  of  both  laminar  and  turbulent  flow  in  a  strongly  curved 
square  duct  with  moderately  thin  boundary  layers  at  the  entrance  to  the  bend.  The 
present  study  concentrates  heavily  on  numerical  solutions  of  the  Navier-Stokes 
equations  for  the  flow  conditions  of  these  measurements  of  Taylor,  Whitelaw  and 
Yianneskis . 

Flow  in  Pipe  Bends 

Rowe  [10]  has  taken  total  pressure  measurements  for  turbulent  flow  in  a  pipe 
bend  of  small  curvature.  For  comparison,  he  also  performed  flow  calculations  based 
on  the  Squire-Winter  inviscid  secondary-flow  approximation.  Patankar,  Pratap  and 
Spalding  [11,  12]  have  performed  calculations  for  both  laminar  and  turbulent  flow 
in  pipe  bends  of  small  curvature,  using  their  "parabolic  flow"  calculation  procedure 
and  a  two-equation/wall-function  turbulence  model.  Agrawal,  Talbot  and  Gong  [13] 
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have  obtained  detailed  measurements  of  laminar  flow  development  in  curved  pipes  with 
uniform  velocity  at  entry.  These  measurements  would  be  useful  for  evaluation  of 
numerical  computations. 


THE  PRESENT  APPROACH 


The  present  objective  is  to  explore  and  evaluate  a  method  for  predicting  turbulent 
flows  in  ducts  and  pipes  of  strong  curvature  by  numerical  solution  of  the  Navier-Stokes 
equations.  A  time-dependent  formulation  is  used  as  an  Iterative  means  of  obtaining 
steady  solutions,  and  the  compressible  form  of  the  equations  are  solved  in  the  low 
Mach  number  regime  (M  &  0.05),  which  closely  approximates  an  incompressible  flow. 

The  governing  equations  are  solved  using  a  consistently-split  linearized  block 
implicit  (LBI)  scheme  developed  by  Briley  and  McDonald  [14,  15].  With  proper  treatment 

t 

of  boundary  conditions,  this  algorithm  provides  rapid  convergence  which  is  not  signifi¬ 
cantly  degraded  by  the  extreme  local  mesh  resolution  which  is  necessary  for  the  near¬ 
wall  sublayer  region  in  turbulent  flows.  The  turbulence  model  used  is  a  one-equation 

model  recently  explored  by  Shamroth  and  Gibeling  [15].  This  model  requires  solution 

2 

of  a  single  equation  governing  turbulence  kinetic  energy  q  ,  in  conjunction  with  an 

algebraically  specified  length  scale.  The  turbulent  effective  viscosity  u  is  obtained 

2  1/2  C 

from  the  Prandtl-Kolomogorov  constitutive  relationship  u  =  £(q  )  .  The  model  also 

\  L 

includes  representation  of  the  influence  of  turbulence  Reynolds  number  on  turbulent 
stress  levels  and  provides  for  resolution  of  the  near-wall  viscous  sublayer  region. 

The  present  effort  concentrates  on  the  geometry  and  flow  conditions  for  which 
detailed 'measurements  have  been  obtained  recently  by  Taylor,  Whitelaw,  and  Yianneskis  [9]. 
This  geometry  is  shown  to  scale  in  Fig.  1  and  consists  of  a  square  duct  with  a  90  degree 
circular  arc  bend  and  with  straight  sections  both  upstream  and  downstream  of  the  bend. 

The  ratio  of  bend  radius  R  to  duct  width  W  is  2.3.  The  measurements  were  taken  for 
Reynolds  numbers  based  on  mean  velocity  and  duct  width  of  790  (laminar  flow)  and 
40,000  (turbulent  flow).  In  each  case,  moderately  thin  shear  layers  (20  -  30%  of 
duct  width)  were  present  at  the  start  of  the  bend.  The  present  study  considers  both 
of  these  measured  flows  and  in  addition  the  corresponding  two-dimensional  channel 
flows  having  the  same  Reynolds  numbers ,  shear  layer  thickness  and  ratio  of 'bend 
radius  to  duct  width.  The  two-dimensional  channel  flows  are  computed  to  establish 
their  flow  structure  for  comparison  with  the  three-dimensional  duct  flow  solutions. 

Test  calculations  are  performed  to  establish  grid  independence  and  to  verify  the 
treatment  of  inflow/outflow  conditions  for  the  two-dimensional  channel  flows. 

The  laminar  and  turbulent  duct  flow  solutions  are  compared  with  the  flow  measurements 
and  with  the  channel  flow  results  to  provide  an  evaluation  of  the  flow  predictions 
obtained  by  the  present  computational  procedure.  Finally,  turbulent  flow  in  a  curved 
pipe  with  radius  ratio  of  24  and  Reynolds  number  of  2.36  x  10^  is  computed  and 
compared  with  the  measurements  of  Rowe  [10]. 
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Coordinate  System 


The  compressible  Navier-Stokes  equations  in  general  orthogonal  coordinates  are 
solved  using  analytical  coordinate  data  for  a  system  of  coordinates  aligned  with  the 
duct  geometry.  The  coordinate  system  is  shown  in  Fig.  2  and  consists  of  an  axial 
coordinate  parallel  to  a  curved  duct  centerline  (which  lies  in  the  Cartesian  x-y 
plane),  and  general  orthogonal  coordinates  x^  in  transverse  planes  normal  to  the 
centerline.  If  the  axial  coordinate  x^  denotes  distance  along  the  centerline  and  if 
K(x^)  =  1/R(x^)  is  the  centerline  curvature,  then  the  metric  scale  factor  h^  for  the 
axial  coordinate  direction  is  given  by  h^  -  1  +  K(x^)  AR^,  x^) ,  where  AR  5  r-R  is 
independent  of  x^.  The  transverse  metric  factors  are  given  by  =  h^  =  1  for 

rectangular  (Cartesian)  cross  sections  and  by  hj  *  1,  h^  «*  x^  for  circular  (polar) 
cross  sections.  The  quantity  AR  is  given  by  AR  =  x^  and  AR  *  X2  cos  x^  for  Cartesian 
and  polar  cross  sections,  respectively.  The  centerline  curvature  K  is  discontinuous 
when  a  straight  segment  of  a  duct  joins  a  constant  radius  segment.  To  remove  the 
associated  coordinate  singularity,  the  flow  geometry  is  smoothed  over  an  axial 
distance  of  one  duct  diameter.  This  is  accomplished  using  a  cubic  polynomial  varia¬ 
tion  in  K  which  matches  both  function  value  and  slope  of  K  at  the  end  points  of  the 
smoothing  region. 


Boundary  and  Initial  Conditions 

The  computational  domain  Is  chosen  to  be  a  region  in  the  immediate  vicinity  of 
the  duct  bend  (cf.  Fig.  1)  embedded  within  a  larger  overall  flow  system  upstream  and 
downstream  of  the  bend.  This  choice  of  computational  domain  requires  inflow  and 
outflow  boundary  conditions  which  adequately  model  the  interface  between  the  computed 
flow  and  the  remainder  of  the  flow  system.  The  inf low/out flow  conditions  used  are 
derived  from  an  assumed  flow  structure  and  are  chosen  to  provide: 

(a)  inflow  with  prescribed  stagnation  pressure  and  temperature  in  an 
inviscid  core  region,  and  with  a  given  axial  velocity  profile  shape  in 
shear  layers,  and 

(b)  outflow  with  a  prescribed  static  pressure  distribution. 

These  boundary  conditions  are  chosen  following  consideration  both  of  a  one-dimensional 
inviscid  characteristics  analysis  and  of  the  physical  process  by  which  many  duct  flows 
are  established.  For  subsonic  flow,  a  characteristics  analysis  of  the  one-dimensional 
inviscid  Euler  equations  indicates  that  two  boundary  conditions  are  required  at  an 


inflow  boundary  and  one  additional  boundary  condition  is  required  at  outflow.  Physically 
a  duct  flow  is  often  established  by  supplying  air  of  a  given  stagnation  pressure  and 
temperature  and  exhausting  the  duct  at  a  given  static  pressure.  The  mass  flux  through 
the  duct  may  then  vary  with  time  until  a  steady  state  is  achieved,  at  which  the  mass 
flux  is  determined  as  a  balance  between  these  inf low/ outflow  quantities  and  viscous 
and  thermal  effects  within  the  duct.  By  choosing  stagnation  pressure  and  temperature 
at  inflow  and  static  pressure  at  outflow  as  the  dominant  boundary  conditions,  the 
present  solution  procedure  allows  both  velocity  and  density  to  vary  with  time,  as  is 
consistent  with  this  physical  process.  As  a  consequence,  pressure  waves  can  be  trans¬ 
mitted  through  the  inflow  boundary  during  the  transient  flow  process  and  are  not 
reflected  back  into  the  computational  domain.  The  reflection  of  pressure  waves  at 
an  inflow  boundary  where  velocity  and  density  are  fixed  in  time  has  often  been  cited 
as  a  cause  of  either  instability  or  slow  convergence  in  other  investigations.  The 
specific  treatment  of  initial  and  boundary  conditions  used  here  is  outlined  below. 

The  initial  and  boundary  conditions  are  devised  from  estimates  of  the  potential 
flow  velocity  Uj.(x^,  x^)  for  the  duct,  a  mean  boundary  layer  thickness  6(x^)  for 

shear  layers  on  transverse  duct  walls,  and  finally  from  an  estimate  of  the  blockage 
correction  factor  B(x1)  for  the  core  flow  velocity  due  to  the  boundary  layer  growth. 

The  potential  flow  velocity  is  approximated  as  uniform  flow  in  straight  segments  of 
the  duct  and  as  proportional  to  r*1  in  curved  segments.  The  constants  (R  -R^)/fcn(R  /Rj) 

and  R  2/2  (R-  R2-R  2)  lead  to  a  unit  mass  flux  for  rectangular  and  circular  cross 
o  o 

sections,  respectively.  Here,  R^  and  Rq  are  the  radii  to  the  inner  and  outpr  walls  of 
the  duct,  and  R  =  (R^  +  Rq)/2  as  in  Fig.  1.  Distributions  of  6(x^)  and  BCx^)  are 
determined  by  recourse  to  a  simple  one-dimensional  momentum  integral  analysis  using  a 
fixed  velocity  profile  shape  and  an  approximate  relationship  between  B  and  mean 
displacement  thickness  6*(x^).  Details  of  this  procedure  are  not  important,  as  the 
results  serve  mainly  as  a  convenient  method  of  setting  approximate  initial  conditions. 
Finally,  a  shear  layer  velocity  profile  shape  f(y/5),  o<fsl  is  chosen  for  each  problem, 
where  y  is  a  parameter  indicative  of  distance  from  a  wall.  The  initial  values  of 
velocity  components  u^,  u2»  u^  are  given  by 


u 


Uj.  B(x1)  f  [y/<5(x)  ]  u2  -  u3  *  0 


(1) 


A  reasonably  accurate  estimate  for  the  pressure  drop  which  will  produce  the  desired 
flow  rate  must  be  made  using  any  convenient  source,  such  as  a  Moody  diagram,  data 
correlations,  momentum  integral  analysis,  or  other  computed  results.  A  smooth  axial 
distribution  of  pressure  which  matches  this  pressure  drop  is  then  assigned  and 
adjusted  to  approximate  local  curvature  of  the  flow  geometry.  This  completes  speci¬ 
fication  of  the  initial  conditions.  It  is  noted  that  although  these  initial  conditions 
do  take  into  account  several  relevant  features  of  the  flow,  the  important  effects  of 
strong  secondary  flows  and  their  distortion  of  the  primary  flow  are  completely 
neglected.  The  resulting  initial  flow  is  thus  a  simple  but  relatively  crude  approxi¬ 
mation  to  the  final  flow  field. 

At  the  inflow  boundary,  a  "two-layer"  boundary  condition  is  devised  such  that 
stagnation  pressure  pQ  is  fixed  in  the  core  flow  region  and  an  axial  velocity  profile 
shape  u^/ue  =  f(y/6)  is  set  within  shear  layers.  Here,  ug  is  the  local  edge  velocity 
which  varies  with  time  and  is  adjusted  after  each  time  step  to  the  value  consistent 


with  p  and  the  local  edge  static  pressure  determined  as  part  of  the  solution.  The 

°  2  2  2  2  2  2 
remaining  inflow  conditions  are  3  U2/3n  =  3  u^/3n  =  0  and  3  c^/Sn  =  g  (x2»  x^) , 

where  n  denotes  the  normal  coordinate  direction  and  c  is  pressure  coefficient.  The 

P  2 

quantity  g  is  computed  from  the  initial  conditions  with  c^  defined  as  l-(BUj)  ,  its 
value  from  the  potential  flow  corrected  for  estimated  blockage.  For  outflow  condi¬ 
tions,  the  static  pressure  is  imposed,  and  second  normal  derivatives  of  each  velocity 


component  are  set  to  zero. 


At  no-slip  walls. 


all  velocity  components  are  set  to  zero. 


and  the  remaining  condition  used  is  3p/3n  =  0,  where  p  is  pressure.  The  condition 
3p/3n  =  0  at  a  no-slip  surface  approximates  the  normal  momentum  equation  to  order  Re 


for  viscous  flow  at  high  Reynolds  number.  Finally,  the  three-dimensional  flow  cases 


are  assumed  to  be  symmetric  about  the  plane  containing  the  curved  duct  centerline. 


and  symmetry  conditions  are  imposed  on  this  boundary. 


Governing  Equations  and  Differencing  Procedures 

The  differencing  procedures  used  are  a  straightforward  adaptation  of  those  used 
by  Briley  and  McDonald  [14]  in  Cartesian  coordinates  for  flow  in  a  straight  duct. 

The  compressible  time-dependent  Navier-Stokes  equations  are  written  in  orthogonal 
coordinates  in  the  form  given  by  Hughes  and  Gaylord  [17].  The  first-derivative 
flux  terms  are  written  in  conservation  form,  and  for  economy  the  stagnation  enthalpy 
is  assumed  constant.  The  definition  of  stagnation  enthalpy  and  the  equation  of  state 
for  a  perfect  gas  can  then  be  used  to  eliminate  pressure  and  temperature  as  dependent 
variables,  and  solution  of  the  energy  equation  is  unnecessary.  The  continuity  and 
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three  momentum  equations  are  solved  with  density  and  the  three  velocity  components 

aligned  with  the  coordinates  as  dependent  variables.  Three-point  central  differences 

were  used  for  spatial  derivatives,  and  second-order  artificial  dissipation  terms  are 

added  as  in  [14]  to  prevent  spatial  oscillations  at  high  cell  Reynolds  number.  This 

treatment  lowers  the  formal  accuracy  to  first  order  but  does  not  seriously  degrade 

accuracy  in  representing  viscous  terms  in  thin  shear  layers.  Analytical  coordinate 

transformations  due  to  Roberts  [18]  were  used  to  redistribute  grid  points  and  thus 

improve  resolution  in  shear  layers.  Derivatives  of  geometric  data  were  determined 

analytically  for  use  in  the  difference  equations. 

The  turbulence  model  used  is  a  one-equation  model  recently  explored  by  Shamroth 

and  Gibeling  [15].  This  model  requires  solution  of  a  single  equation  governing 

2 

turbulence  kinetic  energy  q  ,  in  conjunction  with  an  algebraically  specified  length 

2 

scale  £.  The  equation  governing  the  balance  of  turbulence  kinetic  energy,  q  ,  in 

curvilinear  orthogonal  coordinates  was  derived  from  the  Cartesian  tensor  form  of  this 

equation  given  by  Launder  and  Spalding  [19].  For  the  present  application,  a  very 

simple  three-layer  length  scale  was  constructed  with  the  outer  or  wake  scale  determined 

by  a  one-dimensional  estimate  of  the  growth  of  this  wake  length  scale  from  its  value 

on  inlet.  This  growth  rate  was  essentially  obtained  from  the  Von  Karman  momentum 

integral  equation  assuming  the  wake  length  scale  would  grow  roughly  as  the  boundary 

layer  thickness.  Near  the  nearest  wall  the  length  scale  was  assumed  to  vary  in 

accordance  with  Von  Karman's  linear  relationship  £= icy.  In  the  viscous  sublayer  the 

length  scale  was  damped  by  viscous  effects  according  to  Van  Driest's  formulation. 

Finally,  the  turbulent  effective  viscosity  p  was  obtained  from  the  Prandtl-Kolomogorov 

2  1/2  C 

constitutive  relationship  p^  a  £(q  )  .  The  turbulent  viscosity  was  then  supposed 

isotropic  and  the  stress  tensor  in  the  ensemble  averaged  equations  determined  by 
adding  the  turbulent  viscosity  to  the  kinetic  viscosity.  The  turbulent  kinetic  energy 
near  the  wall  was  damped  out  according  to  the  suggestion  of  Shamroth  and  Gibeling  (Ref.  15) 

Split  LBI  Algorithm 

The  numerical  algorithm  used  is  the  consistently-split  "linearized  block  implicit" 
(LBI)  scheme  developed  by  Briley  and  McDonald  [14,  15]  for  systematic  use  in  solving 
systems  of  nonlinear  parabolic-hyperbolic  partial  differential  equations  (PDE's). 

To  illustrate  the  algorithm,  let 

(<f>n+1-4>n)/At  =  BD(4>n+1)  +  (1-8)  D  (<j>n)  (2) 


«r*  v^-hk» 
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It 


approximate  a  system  of  time-dependent  nonlinear  FDE's  (centered  about  tn+BAt)  for 
the  vector  $  of  dependent  variables ,  where  D  is  a  multidimensional  vector  spatial 
differential  operator,  and  t  is  a  discretized  time  variable  such  that  At*tn+^-tn. 

A  local  time  linearization  (Taylor  expansion  about  ij>n)  is  introduced,  and  this  serves 
to  define  a  linear  differential  operator  L  such  that 


D  (<J>n+1)  =  D  (<pn)  +  Ln  (<|>n+1-<j>n)  +  0  (At2) 


Eq.  (2)  can  thus  be  written  as  the  linear  system 


(I  -  $AtLn)  (<j>n+1-<t>n)  =  AtD  (<J.n) 


without  lowering  the  formal  accuracy. 

The  multidimensional  operator  L  is  divided  into  three  "one-dimensional"  sub-operators 
L^Lj+I^+L^  (associated  here  with  the  three  coordinate  directions),  and  Eq.  (4)  is 
split  as  in  the  scalar  development  of  Douglas  &  Gunn  [20]  and  is  written  as 

(I  -  gAtLj11)  (4>*— 4>n)  -  AtD  (4>n)  (5a) 

(I  -  0AtL2n)  (<( *%")  -  *%n  (5b> 

(I  -  BAtL3n)  (<j)***-4>n)  =  (J>**-(t>n  (5c) 

<|>  1  =  4>  +0(At  )  (5d) 

If  spatial  derivatives  appearing  in  L  are  replaced  by  three-point  difference  formulas, 
then  each  step  in  Eqs.  (5a-c)  can  be  solved  by  a  block-tridiagonal  "inversion". 
Eliminating  the  intermediate  steps  in  Eqs.  (5a-d)  results  in 

(I  -  BAtL^)  (I  -  BAtL?n)  (I  -  0AtL3n)  ($n+1-<j>n)  =  AtD(<|>n)  (6) 

which  approximates  Eq.  (4)  to  order  At^.  Complete  derivations  are  given  by  the 
authors  in  [14,  15],  It  is  noted  that  Warming  and  Beam  [21]  have  suggested  an 
alternate  derivation  for  this  and  other  related  algorithms,  based  on  the  approximate 
factorization  approach  of  Yanenko  and  D'Yakonov  [22],  Their  derivation  is  equivalent 
up  to  Eq.  (4)  and  proceeds  immediately  to  Eq.  (6)  by  observing  that  Eq.  (6)  is  a 
"delta  form"  approximate  factorization  of  Eq.  (4).  Eq.  (6)  is  then  solved  by  intro¬ 
ducing  intermediate  quantities  (<j>  -$n)  and  solving  Eqs.  (5a-c). 


pm** • 


COMPUTED  RESULTS 


Extensive  calculations  were  made  for  the  flow  geometry  in  which  detailed  measure¬ 
ments  were  obtained  by  Taylor,  Whitelaw,  and  Yianneskis  [9].  This  geometry  is  shown 
in  Fig.  1  and  consists  of  a  square  duct  with  a  90  degree  circular-arc-bend  and  with 
straight  sections  both  upstream  and  downstream  of  the  bend.  The  ratio  of  bend 
radius  to  duct  width  is  2.3.  The  measurements  were  taken  for  Reynolds  numbers  based 
on  mean  velocity  and  duct  width  of  790  (laminar  flow)  and  40,000  (turbulent  flow). 

In  each  case,  moderately  thin  shear  layers  (20-30%  of  duct  width)  were  present  at  the 
start  of.  the  bend.  Both  of  these  measured  flows  were  computed.  In  addition  the 
corresponding  two-dimensional  channel  flows  having  the  same  curvature,  Reynolds 
number  and  shear  layer  thickness  were  computed  for  comparison  with  the  three- 
dimensional  duct  flow  solutions. 

Mesh  Refinement  and  Other  Validation  Tests 

Test  calculations  were  first  performed  for  the  two-dimensional  channel  flows 
to  establish  grid  independence  and  to  verify  the  treatment  of  inflow/outflow  conditions. 
The  results  are  shown  in  Figs.  3-5.  In  Fig.  3,  axial  velocity  profiles  at  the  60° 
location  are  shown  for  both  laminar  and  turbulent  flow  and  for  three  radial  grid 
densities.  With  the  26  x  28  grid,  the  mesh  spacing  (in  channel  widths)  adjacent 
to  the  wall  is  0.0026  (turbulent)  and  0.0087  (laminar).  The  results  in  Fig.  3  indicate 
that  26  radial  grid  points  are  sufficient  to  resolve  this  flow.  In  Fig.  4,  the  wall 
pressure  coefficient  c^  is  shown  for  both  laminar  and  turbulent  flow  and  for  two 
axial  grid  densities.  Here  and  elsewhere  in  this  report,  the  axial  coordinate  is 
given  in  degrees  of  turning  within  the  bend  and  in  duct  widths  upstream  or  downstream 
of  the  bend.  The  results  in  Fig.  4  are  not  sensitive  to  the  axial  grid  used.  Note 
that  at  the  inflow  boundary,  a  transverse  pressure  gradient  is  permitted  by  the  imposed 
inflow  boundary  conditions  and  does  occur  in  Fig.  4.  The  static  pressure  is  assumed 
constant  at  the  downstream  outflow  boundary.  In  Fig.  5,  the  effect  of  extending  the 
straight  segments  to  a  length  of  3.5  channel  widths  both  upstream  and  downstream  of 
the  bend  is  shown.  Although  the  extensions  had  no  observable  effect  on  the  turbulent 
flow  results,  there  is  a  noticeable  effect  near  the  90°  location  for  the  laminar  flow. 
This  sensitivity  of  the  laminar  flow  can  be  explained  by  the  occurrence  of  a  very 
small  flow  separation  on  the  inner  wall  in  this  region.  Evidently,  although  the 
solutions  are  insensitive  to  the  axial  grid,  the  downstream  boundary  is  apparently 
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located  too  close  to  this  flow  separation  unless  extensions  are  added.  This  flow 
separation  did  not  occur  in  the  three-dimensional  square  duct  analog  of  this  flow, 
however,  and  the  extensions  were  not  included  in  these  latter  calculations. 

On  balance,  it  is  concluded  that  a  26  x  28  grid,  distributed  to  provide  local 
resolution  for  the  near-wall  region,  was  sufficient  to  achieve  grid  independence  for 
the  two-dimensional  channel  cases.  This  same  grid  was  used  for  the  three-dimensional 
duct  calculations,  and  the  mesh  distribution  for  the  third  coordinate  was  chosen  with 
the  two-dimensional  results  as  a  guide.  Taking  advantage  of  symmetry  about  the  plane 
midway  between  the  endwalls,  a  26  x  28  x  13  grid  was  used,  with  mesh  spacing  adjacent 
to  the  endwall  of  0.0028  (turbulent)  and  0.0091  (laminar)  duct  widths.  Since  the  two- 
and  three-dimensional  flows  have  a  different  structure,  the  resulting  grid  does  not 
guarantee  accuracy  in  three  dimensions.  Nevertheless,  knowledge  of  the  two-dimensional 
accuracy  is  very  useful  in  evaluating  the  three-dimensional  results,  for  which 
extensive  mesh  refinement  results  were  not  feasible. 

Finally,  it  is  necessary  to  establish  and  document  both  the  degree  and  rate  of 

convergence  obtained  in  the  computed  solutions.  The  behavior  of  maximum  change  in 

computed  streamwise  velocity  with  nondimensional  time  is  shown  in  Fig.  6  for  the 

two-dimensional  turbulent  channel  flow.  These  results  are  typical  of  those  obtained 

in  both  two-  and  three-dimensions  and  for  both  laminar  and  turbulent  flow.  Also 

shown  in  Fig.  6  is  the  time  step  number.  Here,  the  particle  residence  time  of  6.75 

represents  the  nondimensional  time  required  for  a  particle  travelling  with  the 

reference  (mean)  velocity  to  traverse  the  duct  centerline  (6.75  duct  widths)  from 

inflow  to  outflow.  In  this  calculation,  the  initial  time  step  was  increased  gradually 

from  0.02  to  0.5  in  the  first  15  steps,  held  constant  at  0.5  for  the  next  15  steps, 

and  thereafter  a  sequence  of  5  time  steps  in  equal  logarithmic  increments  between 

0.01  and  0.5  was  used  cyclically.  The  normalized  maximum  increment  in  velocity  shown 

in  Fig.  6  is  shown  only  for  the  maximum  time  step  of  0.5.  It  should  be  emphasized 

that  Fig.  6  would  not  be  meaningful  if  the  increment  were  not  taken  at  a  large  time 

step,  since  the  temporal  increment  can  be  made  arbitrarily  small  by  taking  a  small 
-A 

time  step.  With  10  as  the  criterion,  convergence  to  the  steady  state  was  usually 

obtained  at  around  80  time  steps  in  both  two-  and  three-dimensional  flow.  Residuals 

in  the  governing  equations  were  also  examined  at  convergence  and  noted  to  be  small  in 

comparison  to  the  other  terms  in  the  governing  equations.  As  a  final  observation, 

-A 

with  10  as  the  criterion,  the  plotted  results  exhibited  little  or  no  discernable 
change  with  further  iteration.  The  behavior  observed  in  Fig.  6  does  represent  rapid 
convergence  for  a  problem  of  this  complexity,  and  this  is  attributed  in  part  to  the 
use  of  an  implicit  algorithm  (including  boundary  conditions)  and  in  part  to  the 
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selection  of  inflow/outflow  boundary  conditions.  The  specification  of  total  pressure 
at  inflow  and  static  pressure  at  outflow  allows  pressure  waves  to  be  transmitted 
through  the  inflow  boundary  during  the  transient  process,  and  this  avoids  the 
instability  or  slow  convergence  often  attributed  to  reflection  of  pressure  waves  and 
an  inflow  boundary  with  fixed  velocity  and  density. 

Three-Dimensional  Curved  Duct  and  Pipe  Flows 

Axial  velocity  profiles  for  turbulent  flow  in  both  the  two-dimensional  channel 
and  in  the  symmetry  plane  x^  =  0  of  the  three-dimensional  duct  are  shown  at  selected 
locations  in  Fig.  7.  Also  shown  are  the  duct  flow  measurements  of  Taylor,  Whitelaw  & 
Yianneskis  [9].  The  two-dimensional  turbulent  channel  flow  was  not  compared  with 
experimental  measurements,  although  as  mentioned  they  were  shown  to  be  essentially 
grid  independent.  The  computed  flow  structure  for  the  channel  flow  consists  of  an 
inviscid  core  flow  with  acceleration  near  the  inner  wall  of  the  bend,  and  with  shear- 
layer  velocity  profiles  characteristic  of  turbulent  flow.  For  comparison,  also  shown 
in  Fig.  7  are  velocity  profiles  in  the  symmetry  plane  for  the  three-dimensional  duct 
flow.  Here,  there  is  evidence  of  considerable  flow  distortion  as  a  result  of  secondary 
flow,  which  is  directed  toward  the  inner  wall  near  the  endwalls  (x^  *  +  0.5)  and  toward 
the  outer  wall  near  the  symmetry  plane  (x^  =  0) .  Although  the  axial  velocity  in  the 
symmetry  plane  initially  accelerates  near  the  inner  wall  up  to  the  30°  location 
(duplicating  the  channel  flow  behavior),  beyond  this  point  the  peak  velocity  shifts 
toward  the  outer  wall.  This  pattern  is  also  present  in  the  data  of  Taylor,  Whitelaw  and 
Yianneskis  [9],  and  the  level  of  agreement  shown  in  Fig.  7  is  generally  good  except 
near  the  inner  wall  near  the  end  of  the  bend.  Radial  velocity  profiles  at  the 
77.5  degree  location  are  shown  in  Fig.  8  and  display  a  very  strong  secondary  flow  near 
the  endwall  with  peak  velocity  up  to  40%  of  the  mean  axial  velocity,  located  very  near 
the  endwall  surface.  Near  the  symmetry  plane,  the  radial  velocity  is  of  order  20%  of 
the  mean  axial  velocity  and  is  directed  toward  the  outer  wall.  The  computed  radial 
velocity  is  in  very  good  agreement  with  the  measurements  except  in  the  region  near 
the  symmetry  plane  and  inner  wall.  On  balance,  the  level  of  agreement  between  predic¬ 
tion  and  experiment  is  good;  the  discrepancy  may  be  due  to  the  turbulence  model,  but 
other  factors  such  as  grid  resolution  may  still  be  present.  In  fact,  for  a  three- 
dimensional  sublayer-resolved  turbulent  flow  calculation  with  relatively  coarse  grid, 
the  level  of  agreement  is  perhaps  better  than  might  be  anticipated. 


To  explore  further  the  accuracy  of  the  predictions,  laminar  flow  in  the  same 
geometry  was  also  considered.  Computed  results  for  laminar  axial  and  radial  velocity 
profiles  corresponding  to  those  given  for  turbulent  flow  are  shown  in  Figs.  9  and  10. 

The  two-dimensional  laminar  flow  in  Fig.  9  resembles  the  corresponding  turbulent  flow 
in  that  the  flow  accelerates  near  the  inner  wall  of  the  bend.  The  velocity  profiles 
are  characteristic  of  laminar  flow,  however,  and  since  the  growth  rate  of  laminar  shear 
layers  is  greater  than  that  of  turbulent  flow,  there  is  less  evidence  of  an  inviscid 
core  region  in  the  downstream  profiles.  In  addition,  laminar  flow  is  less  resistant 
to  separation  in  an  adverse  pressure  gradient,  and  there  is  a  very  small  flow  separa¬ 
tion  at  the  inner  wall  just  downstream  of  the  +0.25  location  (not  actually  shown  in 
Fig.  9a,  but  also  mentioned  earlier  in  the  discussion  of  Fig.  5).  The  three-dimensional 
laminar  duct  flow  profiles  in  Fig.  9b  display  even  more  distortion  downstream  due  to 
secondary  flow  than  the  turbulent  case  in  Fig.  7b.  Although  very  good  qualitative 
agreement  with  the  Taylor,  et  al  measurements  is  obtained,  there  is  quantitative  dis¬ 
agreement,  particularly  near  the  30°  location.  The  radial  velocity  profiles  at  77-1/2° 
are  shown  in  Fig.  10,  and  although  qualitative  agreement  is  obtained,  there  is  quanti¬ 
tative  disagreement  similar  to  that  observed  in  Fig.  8  for  turbulent  flow  but  of 
greater  extent.  Since  the  turbulence  model  has  been  removed  as  a  factor  in  the 
laminar  flow  calculations,  and  since  there  is  no  reason  to  suspect  the  experimental 
data,  the  source  of  disagreement  between  prediction  and  experiment  is  believed  to  be 
mesh  related.  The  two-dimensional  mesh  refinement  tests  appear  to  have  been  inadequate 
as  a  guide  for  the  more  complicated  and  extensive  three-dimensional  flow  structure. 

One  additional  set  of  experimental  comparisons  shown  in  Figs.  11  and  12  provides 
further  clarification  of  these  flow  predictions.  These  two  figures  show  the  axial 
development  along  the  duct  centerline  of  the  axial  and  radial  velocity  components 
for  each  of  the  four  solutions  (2D/3D  and  laminar/turbulent).  In  each  of  the  two- 
dimensional  channel  flows,  the  axial  velocity  (Fig.  11)  undergoes  a  mild  streamwise 
acceleration  due  mainly  to  blockage  caused  by  shear  layer  growth,  and  this  acceleration 
is  larger  in  the  laminar  flow  because  of  the  higher  growth  rate  of  laminar  flow  at 
these  respective  Reynolds  numbers.  In  the  three-dimensional  duct  flows,  the  initial 
(inlet)  velocity  is  higher  than  the  corresponding  channel  flow  because  the  endwall 
shear  layers  of  the  duct  flow  (not  present  in  the  channel  flow)  require  higher 
blockage  to  give  the  same  mean  flowrate.  The  subsequent  downstream  decreases  in 
axial  velocity  for  the  duct  flows  are  a  direct  consequence  of  the  radial  secondary 
flow,  which  convects  low  velocity  fluid  toward  the  centerline.  This  same  behavior  is 
also  present  in  the  experimental  measurements,  although  the  distortion  due  to 
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secondary  flow  is  less  severe  in  the  calculation  than  that  observed  experimentally. 

In  Fig.  12,  the  corresponding  radial  velocities  are  shown  with  a  sign  convention 
such  that  negative  velocity  denotes  flow  toward  the  inner  wall  of  the  bend.  In  the 
two-dimensional  channel  geometry,  the  behavior  is  very  similar  for  both  laminar  and 
turbulent  flow.  There  is  radial  flow  toward  the  inner  wall  near  the  start  of  the 
bend  (0°)  and  toward  the  outer  wall  near  the  bend  exit  (90°).  In  the  three-dimensioanl 
duct  geometry,  the  flow  is  primarily  radially  outward  within  the  bend,  and  the  radial 
flow  is  significantly  stronger  in  the  laminar  case,  as  is  consistent  with  the  develop¬ 
ment  of  the  axial  velocity  in  Fig.  11.  This  same  flow  behavior  is  evident  in  the 
experimental  measurements,  and  although  the  peak  radial  velocity  is  underpredicted 
in  the  laminar  case,  the  agreement  with  experiment  is  remarkably  good  in  the  turbulent 
case.  The  stronger  secondary  flows  observed  in  the  laminar  case  are  believed  to  be 
the  result  of  the  different  axial  velocity  distributions  entering  the  bend.  The 
laminar  velocity  distribution  contains  larger  transverse  (radial)  vorticity  than  the 
turbulent  case  except  in  the  viscous  dominated  region  near  the  walls,  and  based  on 
secondary  flow  theory,  this  should  lead  to  stronger  secondary  flows.  The  laminar 
case  thus  proved  to  be  a  more  difficult  prediction  than  the  turbulent  case. 

Finally,  total  pressure  contours  computed  for  turbulent  flow  in  a  mildly  curved 

circular  pipe  are  compared  with  the  experimental  measurements  of  Rowe  [10]  in  Fig.  9. 

1  2 

Shown  are  contours  of  total  pressure  difference  (pt~pr)  /-^  Prur  *  where  pt  is  total 
pressure  and  p^.,  p^.,  and  ur  are  reference  values  of  static  pressure,  density,  and 
velocity,  respectively,  taken  as  their  values  at  the  pipe  centerline  at  0°.  The 
total  pressure  contours  are  a  sensitive  indicator  of  the  distortion  which  occurs  as  a 
result  of  both  viscous  losses  and  secondary  flow.  In  general  there  is  good  qualita¬ 
tive  agreement  between  prediction  and  measurement.  Again,  the  quantitative  disagree¬ 
ment  which  is  present  is  believed  to  be  mesh-related.  A  30  x  19  x  15  grid  was  used 
for  the  axial,  circumferential  and  radial  directions,  respectively.  The  minimum 
radial  mesh  increment  at  the  pipe  wall  is  0.000765  diameters,  and  the  maximum  radial 
mesh  increment  at  the  centerline  is  0.12  diameters.  This  high  degree  of  mesh 
nonuniformity  was  needed  to  provide  viscous  sublayer  resolution  at  the  Reynolds  number 
(Re.  -  2.36  x  10^)  of  this  case;  the  use  of  only  15  radial  grid  points  under  these 
circumstances  may  be  inadequate. 
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SUMMARY  AMD  CONCLUSIONS 


A  solution  methodology  is  described  for  treatment  of  flow  in  curved  ducts  and 
pipes,  and  it  is  demonstrated  that  rapid  convergence  ( as  80  noniterative  time-steps) 
can  be  obtained  for  three-dimensional  turbulent  flows,  with  sublayer  resoltuion. 

A  series  of  solutions  for  both  laminar  and  turbulent  flow  and  for  both  two  and 
three-dimensional  geometries  of  the  same  curvature  are  presented.  The  accuracy  of 
these  solutions  is  explored  by  mesh  refinement  and  by  comparison  with  experiment. 

In  summary,  good  qualitative  and  reasonable  quantitative  agreement  between 
solution  and  experiment  was  obtained.  Collectively,  this  sequence  of  results  serves 
to  clarify  the  physical  structure  of  these  flows  and  hence  how  grid  selection  procedures 
might  be  adjusted  to  improve  the  numerical  accuracy  and  experimental  agreement.  For  a 
three-dimensional  flow  of  considerable  complexity,  the  relatively  good  agreement  with 
experiment  obtained  for  the  turbulent  flow  case  despite  a  coarse  grid  must  be  regarded 
as  encouraging.  These  results  seem  to  warrant  further  study  to  clarify  the  sources  of 
error  present  and  also  to  perform  calculations  for  additional  geometric  configurations 
and  flow  conditions. 
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(a)  Plane  of  Centerline  (b)  Plane  Normal  to 

Curvature  Centerline 


Fig.  2  -  Schematic  of  Coordinate  System. 
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Axial  Mesh  Refinement  for  Two-Dimensional  Channel 
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Axial  Location 

Fig.  5  -  Effect  of  Extending  Upstream  and  Downstream  Segments  of 

Two-Dimensional  Channel. 
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Radial  Coordinate 


Axial  Velocity,  u 


Transverse  Coordinate 


-  Computed  O  Data  of  Taylor,  et  al 


Radial  Velocity  at  77-1/29 

Fig.  8  -  Radial  Velocity  Profiles  for  Three-Dimensional  Turbulent 

Duct  Flow  (Inner  Wall  to  Right). 


25 


Radial  Coordinate 


(b)  Three-Dimensional  Square  Duct 


Fig.  9  -  Axial  Velocity  Profiles  for  Laminar  Flow  in  Two  and  Three 
Dimensions  (Curves  are  Labeled  by  Axial  Location; 

Xj  *  -0.5  is  Inner  Wall). 
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Radial  Velocity  u2  at  77-1/2° 

-  Radial  Velocity  Profiles  for  Three-Dimensional  Laminar  Duct 
Flow  (Inner  Wall  to  Right). 
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Fig.  11  -  Axial  Development  of  Axial  Velocity  at  Duct  Centerline. 
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Radial  Velocity  u«  at  Duct  Centerline 


-  Total  Pressure  Contours  for  Turbulent  Flow 
in  a  Pipe  Bend  with  Re  -  2.36  x  105,  R/D  -  12 
(Inside  of  Bend  is  to  Right). 
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Dr.  Phillip  S.  Klebanoff 
National  Bureau  of  Standards 
Mechanics  Section 
Washington ,  DC  2023b 

Dr.  G.  Kulin 

National  Bureau  of  Standards 
Mechanics  Section 
Washington ,  DC  2023b 

Dr.  J.  0.  Elliot 

Naval  Research  Laboratory 

Code  8310 

Washington,  DC  20375 

Mr.  R.  J.  Hansen 

Naval  Research  Laboratory 

Code  8bbl 

Washington,  DC  20375 

Professor  A.  Roshko 
California  Institute  of  Technology 
Graduate  Aeronautical  Laboratories 
Pasadena,  CA  91125 

Dr.  Leslie  M.  Mack 
Jet  Propulsion  Laboratory 
California  Institute  of  Technology 
Pasadena,  CA  91103 


Professor  Frederick  X.  Erovend 
University  of  Southern  California 
University  Park 

Department  of  Aerospace  Engineering 
Los  Angeles,  CA  90007 

Professor  John  Laufer 
University  of  Southern  California 
University  Perk 

Department  of  Aerospace  Engineering 
Los  Angeles,  CA  90007 

Professor  T.  P..  Thomas 

Teesside  Polytechnic 

Department  of  Mechanical  Engineering 

Middlesbrough  TS1  3BA,  England 

Dr.  Arthur  B.  Metzner 
University  of  Delaware 
Department  of  Chemical  Engineering. 
Kevark,  DZ  19711 

Professor  Kerry  E.  Bauch 
The  Graduate  School  and  University 
Center  of  the  City  University  of 
Dev  York 

Graduate  Center:  33  West  ^2  Street 
Kev  York,  IIY  10036 

Mr.  IJornan  M.  Nilsen 
Dyntec  Company 

5301  Laurel  Canyon  Blvd.,  Suite  201 
North  Kollyvood,  CA  91607 

Frofessor  L.  Gary  Leal 
California  Institute  of  Technology 
Division  of  Chemistry  and  Chemical 
Engineering 
Fasacena,  CA  91125 

Professor  Robert  E.  Falco 
Michigan  State  University 
Department  of  Mechanical  Engineering 
East  Lansing,  MI  1*SS2L 

Professor  E.  Pune  Lindgren 
University  of  Florida 
Department  of  Engineering  Sciences 
231  Aerospace  Engineering  Building 
Gainesville,  FL  32611 


Technical  Library 
Naval  Missile  Center 
Point  Mugu,  CA  930^1 

Professor  Frencis  R.  Hema 
Princeton  University 
Department  of  Mechanical  and 
Aerospace  Engineering 
Princeton,  NJ  085^0 

Dr.  Joseph  H.  Clarke 
Brovn  University 
Division  of  Engineering 
Providence,  RI  02912 

Professor  J.  T.  C.  Liu 
Brovn  University 
Division  of  Engineering 
Providence,  RI  02912 

Chief,  Document  Section 
Redstone  Scientific  Information  Cent 
Army  Missile  Command 
Redstone  Arsenal,  AL  35809 

Dr.  Jack  V.  Koyt 

Ravel  Ocean  Systems  Center 

Code  2501 

San  Diego,  CA  92152 

Professor  Richard  L.  Pfeffer 
Florida  State  University 
Geophysical  Fluid  Dynamics  Institute 
Tallahassee,  FL  32306 

Dr.  Gary  Chapman 
Ames  Research  Center 
Mail  Step  227-4 
Moffett  Field,  CA  94035 

Dr.  S.  Beus 

Bettis  Atomic  Power  Laboratory 
P.0.  Box  79 

West  Mifflin,  PA  15122 

Dr.  M.  Lubert 
General  Electric  Company 
Knolls  Atomic  Power  Laboratory 
P.0.  Box  1072 
Schnectady,  NY  12301 
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Dr.  A.  K.  M.  Fazle  Hussain 
University  of  Houston 
Department  of  Mechanical  Engineering 
Houston,  TX  7700* 

Professor  John  L.  Lumley 
Cornell  University 
Sibley  School  of  Mechenicel 
and  Aerospace  Engineering 
Ithaca,  KY  1*853 

Professor  K.  E.  Shuler 
University  of  California,  San  Diego 
Department  of  Chemistry 
La  Jolla,  CA  92093 

Dr.  E.  W.  Kontroll 
Fnysical  Dynamics,  Inc. 

P.  0.  Bor.  556 
La  Jolla,  CA  92038 


Professor  Patrick  Leehey 
Massachusetts  Institute  of  Technology 
Department  of  Ocean  Engineering 
Cambridge,  MA  02139 

Professor  Eli  Reshotko 
Case  Western  Reserve  University 
Department  of  Mechanical  and 
Aerospace  Engineering 
Cleveland,  OH  **106 

Professor  P.  S.  Virk 
Massachusetts  Institute  of  Technology 
Department  of  Chemical  Engineering 
Cambri  dge ,  MA  02139 

Professor  E.  Mollo-Christensen 
Massachusetts  Institute  of  Technology 
Department  of  Meteorology 
-/room  5*-1722 
Cambridge,  MA  02139 

Professor  V.  W.  Villmerth 
The  University  of  Michigan 
Department  of  Aerospace  Engineering 
Ann  Arbor,  KI  *8109 

Office  of  Karel  Research 
Code  *8l 

6 00  K.  Quincy  Street 
Arlington,  VA  22217 


Professor  Richard  U.  Miksad 
The  University  of  Texas  at  Austin 
Department  of  Civil  Engineering 
Austin,  TX  78712 

Professor  Stanley  Corrsin 
The  Johns  Hopkins  University 
Department  of  Mechanics  and 
Materials  Sciences 
Baltimore,  MD  21218 

Mr.  M.  Keith  Ellingsvorth 
Power  Program 
Office  of  Naval  Research 
800  N.  Quincy  Street 
Arlington,  VA  22217 

Professor  J.  A.  C.  Humphrey 
Department  of  Mechanical  Engineering 
University  of  California,  Berkeley 
Berkeley,  CA  94720 

Professor  Brian  Launder 

Thermodynamics  and  Fluid  Mechanics  Division 

University  of  Manchester 

Institute  of  Science  &  Technology 

P088  Sackville  Street 

Manchester  M601QD,  England 

Dr.  Simion  Kuo 
Chief,  Energy  Systems 
Energy  Research  Laboratory 
United  Technology  Research  Center 
East  Hartford,  CT  06108 
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